function [finalMatrix] = MatrixPolynomial(p,nmonthvar,lagOrder,nftot,Rho,originalMatrix)

lambdaLagOrder = max(lagOrder+1,p);
numberOfVariables = nmonthvar;

highDMatrix = NaN(size(originalMatrix,1),floor(size(originalMatrix,2)/nftot),nftot);
for j = 1:nftot
    
    % Initialize input matrix
    inputMatrix = originalMatrix(:,j:nftot:end);

    % Initialize shift matrix
    shiftedMatrix = inputMatrix;
    outputMatrix = inputMatrix;

    % Run computation
    for l = 1:lagOrder
        shiftedMatrix = [zeros(numberOfVariables,1) shiftedMatrix(:,1:end-1)];
        repmatMatrix = repmat(Rho(:,l),1,lambdaLagOrder);
        orderMatrix = repmatMatrix.*shiftedMatrix;
        tempMatrix = outputMatrix - orderMatrix;
        outputMatrix = tempMatrix;
    end
    highDMatrix(:,:,j) = outputMatrix;
end

% Compose the matrix back together
finalMatrix = NaN(size(originalMatrix));


for c = 1:size(highDMatrix,2)
    index = nftot*(c-1)+1;
    finalMatrix(:,index:index+nftot-1) = squeeze(highDMatrix(:,c,:));
end




